home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / LinAlg 3.1 / LinAlg / matrix2.cc < prev    next >
Encoding:
Text File  |  1995-12-20  |  6.5 KB  |  198 lines  |  [TEXT/ttxt]
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *            Linear Algebra Package
  6.  *
  7.  *        Basic linear algebra operations, level 2
  8.  *           Matrix transformations and multiplications
  9.  *               of various types
  10.  *
  11.  * $Id: matrix2.cc,v 3.0 1995/01/01 22:29:04 oleg Exp $
  12.  *
  13.  ************************************************************************
  14.  */
  15.  
  16. #include "LinAlg.h"
  17. #include <math.h>
  18. #include <alloca.h>
  19.  
  20. /*
  21.  *------------------------------------------------------------------------
  22.  *                Transpositions
  23.  */
  24.                 // Transpose a matrix
  25. void Matrix::_transpose(const Matrix& prototype)
  26. {
  27.   prototype.is_valid();
  28.   allocate(prototype.ncols,prototype.nrows,
  29.        prototype.col_lwb,prototype.row_lwb);
  30.  
  31.   register REAL * rsp = prototype.elements;    // Row source pointer
  32.   register REAL * tp = elements;
  33.  
  34.                 // (This: target) matrix is traversed in the
  35.                 // natural, column-wise way, whilst the source
  36.                 // (prototype) matrix is scanned row-by-row
  37.   while( tp < elements + nelems )
  38.   {
  39.     register REAL * sp = rsp++;        // sp = @ms[j,i] for i=0
  40.                     // Move tp to the next elem in the col,
  41.     while( sp < prototype.elements + prototype.nelems )
  42.        *tp++ = *sp, sp += prototype.nrows; // sp to the next el in the curr row
  43.   }
  44.   assert( tp == elements + nelems && 
  45.       rsp == prototype.elements + prototype.nrows );
  46.  
  47. }
  48.  
  49. /*
  50.  *------------------------------------------------------------------------
  51.  *            General matrix multiplications
  52.  */
  53.  
  54.             // Compute target = target * source inplace.
  55.             // Strictly speaking, it can't be done inplace,
  56.             // though only the row of the target matrix needs
  57.             // to be saved.
  58.             // "Inplace" multiplication is only possible
  59.             // when the 'source' matrix is square
  60. Matrix& Matrix::operator *= (const Matrix& source)
  61. {
  62.   is_valid();
  63.   source.is_valid();
  64.  
  65.   if( row_lwb != source.col_lwb || ncols != source.nrows ||
  66.       col_lwb != source.col_lwb || ncols != source.ncols )
  67.     info(), source.info(),
  68.     _error("matrices above are unsuitable for the inplace multiplication");
  69.  
  70.                       // One row of the old_target matrix
  71.   REAL * const one_row = (REAL *)alloca(ncols*sizeof(REAL));
  72.   const REAL * one_row_end = &one_row[ncols];
  73.  
  74.   register REAL * trp = elements;    // Pointer to the i-th row
  75.   for(; trp < &elements[nrows]; trp++)    // Go row-by-row in the target
  76.   {
  77.     register REAL *wrp, *orp;               // work row pointers
  78.     for(wrp=trp,orp=one_row; orp < one_row_end;)
  79.       *orp++ = *wrp, wrp += nrows;        // Copy a row of old_target
  80.  
  81.     register REAL *scp=source.elements;        // source column pointer
  82.     for(wrp=trp; wrp < elements+nelems; wrp += nrows)
  83.     {
  84.       register double sum = 0;            // Multiply a row of old_target
  85.       for(orp=one_row; orp < one_row_end;)    // by each col of source
  86.     sum += *orp++ * *scp++;            // to get a row of new_target
  87.       *wrp = sum;
  88.     }
  89.   }
  90.  
  91.   return *this;
  92. }
  93.  
  94.             // General matrix multiplication
  95.             // Create a matrix C such that C = A * B
  96.             // Note, matrix C needs to be allocated
  97. void Matrix::_AmultB(const Matrix& A, const Matrix& B)
  98. {
  99.   A.is_valid();
  100.   B.is_valid();
  101.  
  102.   if( A.ncols != B.nrows || A.col_lwb != B.row_lwb )
  103.     A.info(), B.info(),
  104.     _error("matrices above cannot be multiplied");
  105.  
  106.   allocate(A.nrows,B.ncols,A.row_lwb,B.col_lwb);
  107.  
  108.   register REAL * arp;            // Pointer to the i-th row of A
  109.            REAL * bcp = B.elements;    // Pointer to the j-th col of B
  110.   register REAL * cp = elements;    // C is to be traversed in the natural
  111.   while( cp < elements + nelems )    // order, col-after-col
  112.   {
  113.     for(arp = A.elements; arp < A.elements + A.nrows; )
  114.     {
  115.       register double cij = 0;
  116.       register REAL * bccp = bcp;        // To scan the jth col of B
  117.       while( arp < A.elements + A.nelems )    // Scan the i-th row of A and
  118.     cij += *bccp++ * *arp, arp += A.nrows;    // the j-th col of B
  119.       *cp++ = cij;
  120.       arp -= A.nelems - 1;            // arp points to (i+1)-th row
  121.     }
  122.     bcp += B.nrows;            // We're done with j-th col of both
  123.   }                    // B and C. Set bcp to the (j+1)-th col
  124.  
  125.   assert( cp == elements + nelems && bcp == B.elements + B.nelems );
  126. }
  127.  
  128.             // Compute C = A*B
  129.             // The same is above, only matrix C is already
  130.             // allocated, and it is *this
  131. void Matrix::mult(const Matrix& A, const Matrix& B)
  132. {
  133.   A.is_valid();
  134.   B.is_valid();
  135.   is_valid();
  136.   
  137.   if( A.ncols != B.nrows || A.col_lwb != B.row_lwb )
  138.     A.info(), B.info(),
  139.     _error("matrices above cannot be multiplied");
  140.   if( nrows != A.nrows || ncols != B.ncols ||
  141.       row_lwb != A.row_lwb || col_lwb != B.col_lwb )
  142.     A.info(),B.info(),info(),
  143.     _error("product A*B is incompatible with the given matrix");
  144.  
  145.   register REAL * arp;            // Pointer to the i-th row of A
  146.            REAL * bcp = B.elements;    // Pointer to the j-th col of B
  147.   register REAL * cp = elements;    // C is to be traversed in the natural
  148.   while( cp < elements + nelems )    // order, col-after-col
  149.   {
  150.     for(arp = A.elements; arp < A.elements + A.nrows; )
  151.     {
  152.       register double cij = 0;
  153.       register REAL * bccp = bcp;        // To scan the jth col of B
  154.       while( arp < A.elements + A.nelems )    // Scan the i-th row of A and
  155.     cij += *bccp++ * *arp, arp += A.nrows;    // the j-th col of B
  156.       *cp++ = cij;
  157.       arp -= A.nelems - 1;            // arp points to (i+1)-th row
  158.     }
  159.     bcp += B.nrows;            // We're done with j-th col of both
  160.   }                    // B and C. Set bcp to the (j+1)-th col
  161.  
  162.   assert( cp == elements + nelems && bcp == B.elements + B.nelems );
  163. }
  164.  
  165.             // Create a matrix C such that C = A' * B
  166.             // In other words,
  167.             // c[i,j] = SUM{ a[k,i] * b[k,j] }
  168.             // Note, matrix C needs to be allocated
  169. void Matrix::_AtmultB(const Matrix& A, const Matrix& B)
  170. {
  171.   A.is_valid();
  172.   B.is_valid();
  173.  
  174.   if( A.nrows != B.nrows || A.row_lwb != B.row_lwb )
  175.     A.info(), B.info(),
  176.     _error("matrices above are unsuitable for A'B multiplication");
  177.  
  178.   allocate(A.ncols,B.ncols,A.col_lwb,B.col_lwb);
  179.  
  180.   register REAL * acp;            // Pointer to the i-th col of A
  181.            REAL * bcp = B.elements;    // Pointer to the j-th col of B
  182.   register REAL * cp = elements;    // C is to be traversed in the natural
  183.   while( cp < elements + nelems )    // order, col-after-col
  184.   {
  185.     for(acp = A.elements; acp<A.elements + A.nelems;)
  186.     {                    // Scan all cols of A
  187.       register double cij = 0;            
  188.       register REAL * bccp = bcp;        // To scan the jth col of B
  189.       for(register int i=0; i<A.nrows; i++)    // Scan the i-th row of A and
  190.     cij += *bccp++ * *acp++;            // the j-th col of B
  191.       *cp++ = cij;
  192.     }
  193.     bcp += B.nrows;            // We're done with j-th col of both
  194.   }                    // B and C. Set bcp to the (j+1)-th col
  195.  
  196.   assert( cp == elements + nelems && bcp == B.elements + B.nelems );
  197. }
  198.